Municipal Crime in Residential Subdivisions

\(\color{red}{\text{Notes for Group 6}}\)

0.1 Group and Final Project Deadlines

• Reports Final Drafts: \(\color{red}{\text{May 16, 2021 Sunday 11:00pm EST}}\)
• Report Merge: May 17, 2021 Monday EST
• Youtube Recording Prepared and Uploaded: May 19, 2021 Wed
• Project Due Date: \(\color{red}{\text{Submitted by May 20, 2021 Thursday 11:59pm EST}}\)

Please don’t run tuning code or calibration code in the final version of your result. Just tune the model and give the final working version for inclusion. You can either comment out or add an eval=F flag to the relevant code chunk for such code.

0.2 Components

The raw data files are used to generate the cleaned data files. The cleaned data files are: crimesdiv.geojson and div.geojson. You will be able to use them for your data science work and ignore the 4 other raw files.

1. Exploratory Data Analysis \(\color{red}{\text{Randy}}\)

Visualizations are needed

• Overlay map of census tracts and Town of Cary boundary limits with background tile from real map
• Map of weighted average income by residential subdivision with annotations/title.
• Map of total number of crimes by residential subdivision with annotations/title.

The eBook Geocomputation with R has a lot of good examples using tmap working with sf. (https://geocompr.robinlovelace.net/adv-map.html#interactive-maps)

• Bar chart - Total Crime count by month from Jan 2018 - Nov 2020
• Heat map - Crime Time by month from Jan 2018 - Nov 2020

Examples of traditional charts and heatmaps specific to crime data (https://rpubs.com/Rjacks14/633448)

• Small multiples of crime maps by subdivision showing different crime categories.
• If possible (optional) Interactive map showing income and crime by subdivision in shiny or leaflet using tmap.

2. Unsupervised Analysis \(\color{red}{\text{Phil}}\)

   k-means clustering analysis to identify the centroid of each region. I suggest mapping each region to its equivalent centroid.

Interesting use of k-means clustering along with data visualization. But the distance metric is badly implemented and I wouldn’t follow this metric but the use of only 2 clusters is interesting.

(https://towardsdatascience.com/exploring-clustering-and-mapping-torontos-crimes-96336efe490f)

I suggest working hard to define a distance metric that includes:

• the type of crime but may bucket the ucr codes into a small number of categories.
• median household income
• try using 2 (high vs low) and (4 districts) and more as needed.
• i.e plain Euclidean distance of the actual points does not seem great.

Cluster the data in 2018 alone. Re-run on 2019 data and see if the cluster centroids are the still the same. Plot both diagrams. Color code the residential subdivisions in a map showing the k-means.

The benefit of k-means is not having to explain statistical significance or why you did not use spatial correlations.

3. Principal Components Analysis \(\color{red}{\text{Alex}}\)

Build a matrix of size 713 X 30 of subdivisions and attributes like crime_type and income. Run a PCA and interprete the first and second components a la UC-R article for US Arrests. (https://uc-r.github.io/pca) Will focus on using 2018 data.

4. Supervised learning: Spatial Regression Model \(\color{red}{\text{Alex & TBD}}\)

Not sure if we can finish this component by target date. The challenge is that regression models require handling spatial auto-correlation. A lot to enjoy!

Rspatial website has a good treatment spatial analysis and the supporting software. (https://rspatial.org/raster/analysis/1-introduction.html)

Bivand, Pebesma, Gomez-Rubio Applied Spatial Data Analysis with R (http://gis.humboldt.edu/OLM/r/Spatial%20Analysis%20With%20R.pdf)

Lansley, Cheshire An Introduction to Spatial Data Analysis and Visualisation in R (https://spatialanalysisonline.com/An%20Introduction%20to%20Spatial%20Data%20Analysis%20in%20R.pdf)

Introduction

This document analyzes municipal crime data for the Town of Cary, North Carolina using machine learning methods. In particular, we examine patterns in crime incidents over crime category and residential subdivision. We use both unsupervised machine learning methods to identify clusters and patterns of crime activity and supervised methods to predict patterns of activities.

You may wish to turn off the code in section 1

Section 1 describes the background and data sources.

Section 2 conducts data wrangling and interpolate data to obtain the dataset for our analysis. We construct and describe datasets in later sections.

Section 3 applies a k-means clustering algorithm to identify crime hot spots

Section 4 examines principal component analysis of the crime data

Section 5 examines spatial regression models to predict crime rates using the 2018 training and 2019 test split.

Section 6 concludes our remarks.

Section 7 presents our R code and technical appendices and references.

1 Background and Data Sources

1.1 Background

Cities and towns have good and bad neighborhoods. One common measure of neighborhood health is the level of crime. Another is household income. We wanted to investigate crime and income inequality by neighborhood in a mid-sized town from a statistical learning approach. Such towns have been less studied than larger metropolitan areas such as San Francisco, Chicago or New York. Such an investigation needs detailed crime and income data and a suitable definition of neighborhoods. It also requires the application of statistical learning methods to geospatial data along with the above mentioned attributes.

We investigate the Town of Cary in Wake County, North Carolina. This is a relatively prosperous middle sized boom town located in the Research Triangle area of North Carolina where technology and health companies are concentrated. The Town of Cary has an open data initiative. On the Town’s website, the Police Department publishes crime incident data compiled by its Police Department, residential subdivision boundary data and town boundaries. Lastly, recent US Census Bureau American Community Survey (ACS) data on household income is available to measure economic affluence.

\(\color{red}{\text{More text to explain structure of the paper goes here.}}\)

1.2 Data Sources

The data sources consists of 4 data files containing typical numerical and qualitative data and geospatial data.

• Town of Cary Boundary File: cary_boundary.geojson which contains the geographic border information of the town.

• Census Income Data: We used the American Community Survey (ACS) dataset for 2018 median household income on census tracts overlapping with the Town of Cary.

• Residential Subdivisions: We use the file: httpmapstownofcary0.geojson obtained from the Town of Cary Open data portal to identify approximately 700 residential subdivisions within the Town. These correspond to neighborhoods of the Town.

• Crime Incident Data from the Cary Police Department: We use the historical crime from the Town of Cary Police Department as of Nov 30, 2020. This file contains all historical crime incidents with geospatial data, date, crime category information. We will use the crime incident data from 2018 to 2020.

First we load the relatively small files containing town and residential subdivision and household income data.

Then we load the crime incident data. The raw dataset is quite large and requires some effort to import.

After importing all the raw files, we can examine their dimensions and contents below:

Raw data file dimensions

	town_limits	acs_income	subdivisions	crimes
number rows	1	47	717	101405
number columns	1	5	8	37
source	NA	NA	NA	NA

1.3 Data Issues and How We Address Them

To conduct this analysis, we had to address both the usual missing or bad data issues but also geospatial and methodological issues. We review each data source and discuss both types of issues jointly.

1.3.1 Town Boundaries

The Town of Cary straddles two counties and a significant part of the Town is zoned for commercial or business purposes. However, we determined that the residential subdivisions lie inside the Town thus avoiding the complex issue where subdivisions are split across two municipalities. This dataset provided in GeoJson format is mainly a control during exploratory data analysis. We rely on the residential subdivisions rather than the town boundary datafile to conduct our statistical learning analyses.

1.3.2 Residential Subdivisions

Residents in Cary typically live in residential subdivisions. These subdivisions are overseen by homeowner associations (HOA) but are not municipal administrative units. Such subdivisions are often built on farmland or previously vacant land lacking basic utilities like water, sewage or power prior to the development. Rather residential subdivisions are planned communities constructed by real estate development firms to be sold to homeowners.

We found and removed 4 subdivisions out of 717 where geometry data is missing. We added an identity column div_id to serve as a primary key because subdivision names did not appear to be unique.

## [1] "raw subdivisions data: 717 vs. tidied subdivisions data: 713"

1.3.3 Census Tract Level Household Income

Next, we obtained all census tracts covering the Town of Cary limits. It was necessary to construct a dataset with 47 census tracts to cover the Town of Cary limits. However, census tracts don’t generally align with residential subdivisions. Therefore, multiple census tracts may partition a residential subdivision. The data was obtained from the non-profit censusreporter.org website as a Geojson file. There was no missing data but the median household income data element was renamed from the obscure B19013001 to income for clarity. The data comes from the 2018 ACS dataset. We display some representative rows below while simplifying the dataset further.

Lastly, the table below shows the median household income in dollars within each census tract based on the ACS 2018 Survey.

geoid	tract_name	income	geometry
14000US37037020701	Census Tract 207.01, Chatham, NC	82833	MULTIPOLYGON (((1991620 760…
14000US37037020702	Census Tract 207.02, Chatham, NC	70938	MULTIPOLYGON (((1959672 696…
14000US37183053003	Census Tract 530.03, Wake, NC	68050	MULTIPOLYGON (((2070867 719…
14000US37183053004	Census Tract 530.04, Wake, NC	102143	MULTIPOLYGON (((2065252 723…
14000US37183053005	Census Tract 530.05, Wake, NC	155673	MULTIPOLYGON (((2060975 711…

Because our goal is to estimate median household income of a residential subdivision using census tract data, the non-alignment of these two geographical datasets poses a challenge. We estimate the median household income \(f(S)\) for a residential subdivision \(S\) by taking a weighted average of the median household income of all census tracts \(c\) that intersect a given residential subdivision. More formally, we define:

\[ f(S) = \sum_{ c \cap S \neq \emptyset } \frac{\mu(c \cap S) }{\mu(S)} \text{Income}(c) \] where \(\mu(A)\) is the area of a region \(A\).

To construct this weighted average requires geospatial intersections and dataframe joins which we do using the sf geospatial library and tidyverse for typical dataframe operations.

We observe that about 6% of all subdivisions 43 out of 713 belong to more than one census tract. Of those, only 4 belong to 3 census tracts. Moreover, the median household income range across the census tracts is not extreme.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   48314   82946  109628  104216  123613  169600       1

The lowest income is $48,314 while the highest is $169,600 with a median of $109,628.

## [1] 756  12

We compute the weighted average income below and then save its value back to the residential subdivision dataset. Sample rows are illustrated below.

div_id	category	name	shape_stlength	shape_starea	approvedlots	description	geo_point_2d	avg_income	geometry
1	Single-Family Detached	Southwick	4109.284	801167.3	49	Single-Family Detached	35.77878, -78.84472	131875.00	MULTIPOLYGON (((2045503 737…
2	Single-Family Detached	Forest Creek	7206.395	1847059.2	92	Single-Family Detached	35.71299, -78.78806	155673.00	MULTIPOLYGON (((2063243 713…
3	Single-Family Detached	Stonehaven	5661.140	805679.6	52	Single-Family Detached	35.75736, -78.76810	95867.19	MULTIPOLYGON (((2067795 730…

1.3.4 Crime Incidents

Crime incident data is a large historical dataset. The raw data file, obtained in Dec 2020, contains crime incident data from 1977 to 2020. We have decided to use only the most recent years (2018-2020).

This time limitation still produces a large number of observations (over 13000) but avoids two problems. One is the data quality issues arising from inconsistent reporting especially during its early years. The second, equally important, is that the Town has grown rapidly from the early 2000s. New residential subdivisions have sprouted up over the last two decades. Thus, the absence of crime in a residential subdivision could mean that it was not yet built. Unfortunately, we don’t have historical data on the creation of residential subdivisions. So we are restricted to doing a point-in-time study.

We also filter the incidents by data quality. A significant number of crime incidents lack geometry data which means we cannot find the location within a subdivision. We purge those observations but they are less than half. In recent years, we find geospatial data to be quite completely.

Lastly, we restrict the crime dataset to those incidents occuring in the residential zones of the town. Roughly forty percent of crime appears to occur in commercial, business or non-residential zones. These are incidents within the town limits but not within any residential subdivision. These are important for quality of life issues within the town but less important for differentiating between residential subdivisions.

## [1] 13732    12

	incident_number	newdate	ucr	crime_type	crime_category	geocode	district	geometry
1	18001463	2018-02-17	9914	ALL OTHER - ALL TRAFFIC EXCEPT DWI (NON-UCR)	ALL OTHER	HAMPTON WOODS LN	D1	POINT (2074540 743896.3)
3000	20001132	2020-02-03	23C	LARCENY - SHOPLIFTING	LARCENY	CROSSROADS BLVD	D3	POINT (2077415 731775)
4000	20002505	2020-03-17	90Z	ORDINANCE - ANIMAL BITE	ALL OTHER	CARPENTER FIRE STATION RD	D2	POINT (2030835 754235)
6000	18000128	2018-01-04	220	BURGLARY - FORCIBLE ENTRY	BURGLARY	PETTY FARM RD	D2	POINT (2036305 758545)

In the data table above, we illustrate some representative rows and columns from the crime data table. We briefly comment on the columns:

• The incident_number serves as the unique key in the crime dataset.
  
• The crime is recorded on newdate.
  
• The ucr code follows the FBI Uniform Crime Reporting Program standards in classifying incidents but also uses a small number of non-standard codes which appear to be unique to the Town.
• The crime_type is the most granular but has a many-to-one relationship to the ucr code.
• The district represents one of the four council districts which partitions the Town.

Finally, we inner join the crime incidents data to the residential subdivision and median household income data previously discussed to obtain our cleansed dataset crimes_div. We can see below that the total number of crimes in 2018 , 2019 and a partial 2020 are quite stable below.

## [1] 6239   21

Finally, we export the enriched crimes data and residential subdivision dataset for re-use later.

1.3.5 Coordinate Transformations

We change the coordinate reference system from a geometric crs to a projected CRS on all datasets. We see EPSG=2264 which is the NAD83/North Carolina coordinate system commonly used by the state agencies. The transition from a geographic CRS to a projected CRS means distances and areas are easier to measure. EPSGO.io (https://epsg.io/2264-15835) states the projection accuracy is within 1.0 meter on a bounding box around the entire state.

1.3.6 Cleansed Crime Data

In this final section, we reload the cleansed data files representing the two dataframes for our modeling work. We load crimesdiv which contains the crime incident data enriched with geographical location, associated residential subdivision and weighted average median household income. We also load the dataset of residential subdivisions div.

## Reading layer `crimesdiv' from data source `/Users/philiptanofsky/Documents/School/CUNY/MSDS/Courses/DATA605/DATA622_MACHINELEARNING/FINAL_PROJ/work/crimesdiv.geojson' using driver `GeoJSON'
## Simple feature collection with 6239 features and 19 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2017383 ymin: 692785 xmax: 2080146 ymax: 769015
## CRS:           2264

## Reading layer `div' from data source `/Users/philiptanofsky/Documents/School/CUNY/MSDS/Courses/DATA605/DATA622_MACHINELEARNING/FINAL_PROJ/work/div.geojson' using driver `GeoJSON'
## Simple feature collection with 713 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2016873 ymin: 691265.6 xmax: 2080271 ymax: 769129.4
## CRS:           2264

2 Visualization and Statistical Summary

3 k-Means Clustering

4 Unsupervised learning

To evaluate the relationship between income and crime rates in Cary, this analysis applies the k-means clustering approach to the given data set. K-means clustering is an unsupervised machine learning algorithm that separates the data into k clusters in which the value of k is supplied beforehand. Each cluster is defined by the centroid, or numeric center of the cluster based on the mean of points comprising the cluster.

The analysis uses the standard algorithm, the Hartigan-Wong algorithm, which calculates the within-cluster variation as the sum of squared distances between the points and the assigned cluster’s centroid. The analysis relies on Euclidean distance to maintain simplicity. For evaluating the k-means model, the total within-cluster sum of squares measures, referred to as the the compactness or goodness of the model, is to be as small as possible.

Initial attempts to use the geospatial data in the k-means algorithm itself proved less than optimal as the geospatial appeared to skew the results which rely on a Euclidean distance metric. Instead, the geospatial data is not used in the k-means model clustering, but later used to plot the crime instances on map based on cluster.

Additionally, early attempts to derive the k-means model based on more than two independent variables proved difficult to interpret as the k-means algorithm uses dimensions to define the clusters. To overcome this challenge of interpretability, the decision was made to only use two independent variables per model in order to provide clear understanding of the relationship between the two variables and clear comprehension of the resulting k clusters. For each model, the subdivision median income and a defined crime count comprise the inputs.

First, the analysis attempts to describe a relationship between the subdivisions of Cary by income and total crime count for the years 2018 and 2019. Then, the cluster models are built off median income and count of assaults or theft by subdivision in attempt to better understand the areas with propensity for crime.

## Reading layer `crimesdiv' from data source `/Users/philiptanofsky/Documents/School/CUNY/MSDS/Courses/DATA605/DATA622_MACHINELEARNING/FINAL_PROJ/work/crimesdiv.geojson' using driver `GeoJSON'
## Simple feature collection with 6239 features and 19 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2017383 ymin: 692785 xmax: 2080146 ymax: 769015
## CRS:           2264

4.1 Data evaluation

For the calendar year 2018, the city of Cary has 2,188 reported crimes occurring in 485 unique subdivisions as defined by the \(div_id\).

## [1] 485

## [1] 2188

Based on data inspection, the incident number is not unique by crime instance. To remove any redundancy of crime counts, the resulting number of unique crimes per incident number is 1,965 for Cary in 2018.

## [1] 1965

Despite the lack of readability due to the large number of subdivisions, the distance matrix below based on the classical method of Euclidean distance measure finds plenty of large dissimilarities as evidenced by red along with fairly similar subdivisions as seen in blue. The distance matrix confirms enough dissimilarity exists to continue with the k-means clustering approach.

4.2 Baseline Model

As an initial attempt of the k-means algorithm, a model of just two clusters is attempted to define a baseline for later models. Note, as the k-means algorithm requires all numeric data, the data is scaled before included in the model to ensure outliers do not weigh too heavily on the model’s outcome as the model is based on Euclidean distance measurements.

## List of 9
##  $ cluster     : Named int [1:485] 2 1 1 2 2 1 1 2 2 1 ...
##   ..- attr(*, "names")= chr [1:485] "2" "4" "8" "9" ...
##  $ centers     : num [1:2, 1:2] -0.925 0.638 0.457 -0.315
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "1" "2"
##   .. ..$ : chr [1:2] "avg_inc" "n"
##  $ totss       : num 968
##  $ withinss    : num [1:2] 404 209
##  $ tot.withinss: num 612
##  $ betweenss   : num 356
##  $ size        : int [1:2] 198 287
##  $ iter        : int 1
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

With a simple baseline model of two clusters, the model results in the following composition.

• cluster: The vector of 485 integers indicating the cluster for each point in the initial dataframe.

• centers: The matrix of the two centers.

• totss: The total sum of squares result is 968.

• withinss: The vector of within-cluster sum of squares results in 209 and 404 for the two clusters.

• tot.withinss: The key metric, the total within-cluster sum of squares is 612.

• betweenss: The between-cluster sum of squares (totss- tot.withinss) is 356.

• size: The number of points in each cluster is 287 and 198.

The standard pairwise scatterplot illustrates the two clusters based on the original variables of average income (\(avg_inc\)) and total crime count (\(n\)). As seen below, the scatterplot does not separate in way to provide a clear interpretation as the cluster do not directly overlap but appear to separate based primarily on average income.

Given the lack of clarity with the two-cluster baseline model, k-means models are generated with cluster sizes of 3-5.

With more clusters, a clear differentiation between each cluster appears improving the ability to define each cluster by average income and total count of crime,

4.3 Selecting k clusters

With the requirement of pre-determined cluster counts before running the model, an assessment of the cluster sizes is needed to select the final k cluster number.

First, the elbow method is used by calculating for each k, 1-10, the total within-cluster sum of square (wss). A clear “elbow” appears at 3 k clusters.

Next, the average silhouette method which suggests an optimal number of 10 k clusters.

Finally, the gap statistic method recommends a size of 1 cluster. Given the goal of creating clusters or categories, 1 is not a reasonable suggestion.

In an attempt to compare k-means algorithms, the PAM clustering function was considered as part of the analysis. The resulting plot of 4 clusters very nearly matches the kmeans clustering function from above. An attempt at 5 clusters with the PAM algorithm only divided the sections along the x axis, in essence using income more than crime count. After comparison, the PAM algorithm did not prove better than kmeans, thus the analysis continues to rely on the kmeans implementation.

Using the 5-cluster model with the kmeans algorithm from above, the \(custplot\) function plots the resulting clusters along with the explanation that the two components explain 100% of the point variability. The 100% result is expected based on the use of only two input variables.

After consideration of the model results and goals of this analysis, the 5-cluster kmeans model is selected despite the suggestions of the optimal cluster evaluations. The 5 clusters can easily be interpreted as:

• Low income, low crime rate

• Medium income, low crime rate

• High income, low crime rate

• Medium crime rate

• High crime rate

## List of 9
##  $ cluster     : Named int [1:485] 4 5 1 5 5 5 1 4 5 2 ...
##   ..- attr(*, "names")= chr [1:485] "2" "4" "8" "9" ...
##  $ centers     : num [1:5, 1:2] -1.153 -0.545 -1.401 1.364 0.124 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:2] "avg_inc" "n"
##  $ totss       : num 968
##  $ withinss    : num [1:5] 34.8 58.4 12.5 52 53.3
##  $ tot.withinss: num 211
##  $ betweenss   : num 757
##  $ size        : int [1:5] 119 53 6 110 197
##  $ iter        : int 4
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

For the final 5-cluster kmeans model, the model definition includes:

• totss: 968 - same as the 2-cluster model as expected.

• withinss: All 5 values below 60, much better than the 209 and 404 from the baseline.

• tot.withinss: 211, a marked improvement from the 612 of the baseline, as this is the key metric for assessing the model.

• betweenss: 757, a result much higher than compared to the baseline model, another improved metric

• size: The five clusters result in counts of 119, 53, 6, 110, 197. Three clusters contain more than 100 instances, one cluster with roughly 50 and a final a outlier cluster with just 6 instances. Despite the incongruity of the cluster counts, the outlier cluster represents the high crime count and thus will help provide clarity in finding the crime hot spots.

4.4 Final Models

Based on the selection of 5 k clusters, the models are generated for the years 2018 and 2019 to compare the year-to-year relationship. With each subdivision labeled for one of the 5 clusters, the initial list of crimes are then mapped to a single cluster based on the subdivision of the crime. The resulting maps of crime occurrences in Cary for 2018 and 2019 portray each crime as color coded for the cluster of the corresponding subdivision.

To better understand each cluster, the centroid values are calculated against the prescaled data. For the year 2018, the results are as follows.

• Cluster 1: Low income ($65950), low crime count (2.79)

• Cluster 2: Medium crime count (11.9)

• Cluster 3: High crime count (28.2)

• Cluster 4: High income ($146763), low crime count (2.92)

• Cluster 5: Medium income ($106938), low crime count (2.60)

For the year 2019, the results are as follows.

• Cluster 1: Low income ($65134), low crime count (2.62)

• Cluster 2: High crime count (21.6)

• Cluster 3: Medium income ($106687), low crime count (2.34)

• Cluster 4: Medium crime count (9.53)

• Cluster 5: High income ($147120), low crime count (2.54)

Both 5-cluster models result in the same five understood clusters, yet the kmeans model does not result in the same order between the two models - an interesting artifact of the kmeans algorithm as each model is generated independently and the results are quite similar.

The below map plots every reported crime in Cary for 2018 with the color corresponding to the category of the subdivision using median income and total crime count as input to the clustering. The map shows the clear distinctions of the low, medium and high income sections of the city. The crime hot spots do appear in the northeastern section of the city. The medium crime count areas do appear interspersed throughout the city. Given a medium median crime count of less than 12, then even in these subdivisions, the count is less than one crime per month.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

The 2019 crime map of Cary based on subdivision category appears quite similar to the 2018 map. Given the difference of one year, the similarity is expected.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

4.5 Specific crimes

In an attempt to tease out the locations for certain types of crimes, the above approach of defining each crime by the corresponding subdivision’s median income and total crime count is applied separately to crimes labeled as assault and theft. The category of assault is comprised of aggravated assault and simple assault. The category of theft is defined by the sum of larceny, stolen property and motor vehicle theft. For the specific crime models, the years 2018 through 2020 are considered.

4.5.1 Evaluation of Assault Counts

Total assaults in Cary by year:

• 2018: 485

• 2019: 474

• 2020: 457

## [1] 485   2

## [1] 474   2

## [1] 457   2

To more easily assess the clustering by count, the tidymodels approach to kmeans clustering generates the plots for k clusters of 1 through 9.

Based on the “elbow method”, 3 clusters is the recommendation.

## List of 9
##  $ cluster     : Named int [1:485] 1 3 3 3 1 3 3 1 3 3 ...
##   ..- attr(*, "names")= chr [1:485] "2" "4" "8" "9" ...
##  $ centers     : num [1:3, 1:2] 0.857 -0.617 -0.775 -0.265 3.402 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Assault"
##  $ totss       : num 968
##  $ withinss    : num [1:3] 123.9 69.2 150.6
##  $ tot.withinss: num 344
##  $ betweenss   : num 624
##  $ size        : int [1:3] 228 25 232
##  $ iter        : int 2
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

## List of 9
##  $ cluster     : Named int [1:474] 1 3 3 3 3 1 3 1 3 3 ...
##   ..- attr(*, "names")= chr [1:474] "2" "3" "4" "8" ...
##  $ centers     : num [1:3, 1:2] 0.905 -0.458 -0.742 -0.27 2.449 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Assault"
##  $ totss       : num 946
##  $ withinss    : num [1:3] 107 142 105
##  $ tot.withinss: num 354
##  $ betweenss   : num 592
##  $ size        : int [1:3] 206 44 224
##  $ iter        : int 3
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

## List of 9
##  $ cluster     : Named int [1:457] 1 3 3 3 3 1 3 3 1 3 ...
##   ..- attr(*, "names")= chr [1:457] "2" "4" "8" "9" ...
##  $ centers     : num [1:3, 1:2] 0.906 -0.771 -0.73 -0.231 2.519 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Assault"
##  $ totss       : num 912
##  $ withinss    : num [1:3] 123.8 83.4 90
##  $ tot.withinss: num 297
##  $ betweenss   : num 615
##  $ size        : int [1:3] 205 44 208
##  $ iter        : int 3
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

Based on the model definition outputs, the resulting tot.withinss for the 3 years:

• 2018: 344

• 2019: 354

• 2020: 297

The results aren’t quite as good as the 5-cluster models of the previous section but to be expected given the lower count of k clusters. The three values are close enough in range to consider the model comparisons valid.

For 2018, the resulting cluster definitions for crime of assault:

• Cluster 1: High income, low assault count

• Cluster 2: High assault count

• Cluster 3: Low income, low assault count

For 2019, the resulting cluster definitions for crime of assault:

• Cluster 1: High income, low assault count

• Cluster 2: High assault count

• Cluster 3: Low income, low assault count

For 2020, the resulting cluster definitions for crime of assault:

• Cluster 1: High income, low assault count

• Cluster 2: High assault count

• Cluster 3: Low income, low assault count

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

Given the small range of three years, the maps for each year appear very similar in identifying the crime hot spots based on assault.

4.5.2 Evaluation of Total Theft Counts

The same approach as above is followed to cluster the theft counts by subdivision in relationship to median income.

Again, the tidymodels approach is used to cluster for sizes of 1-9. As above, the plot for three clusters follows the same visual distinction as the plot for assault counts. Also, the elbow method recommends three clusters.

## List of 9
##  $ cluster     : Named int [1:485] 3 1 1 1 3 1 1 3 1 1 ...
##   ..- attr(*, "names")= chr [1:485] "2" "4" "8" "9" ...
##  $ centers     : num [1:3, 1:2] -0.783 -0.471 0.865 -0.171 2.766 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Theft"
##  $ totss       : num 968
##  $ withinss    : num [1:3] 122 122 130
##  $ tot.withinss: num 374
##  $ betweenss   : num 594
##  $ size        : int [1:3] 227 34 224
##  $ iter        : int 2
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

## List of 9
##  $ cluster     : Named int [1:474] 3 2 2 1 2 1 2 3 2 2 ...
##   ..- attr(*, "names")= chr [1:474] "2" "3" "4" "8" ...
##  $ centers     : num [1:3, 1:2] -0.503 -0.76 0.916 2.216 -0.286 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Theft"
##  $ totss       : num 946
##  $ withinss    : num [1:3] 128.6 99.1 120.1
##  $ tot.withinss: num 348
##  $ betweenss   : num 598
##  $ size        : int [1:3] 52 215 207
##  $ iter        : int 3
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

## List of 9
##  $ cluster     : Named int [1:457] 1 2 3 2 2 1 2 2 1 2 ...
##   ..- attr(*, "names")= chr [1:457] "2" "4" "8" "9" ...
##  $ centers     : num [1:3, 1:2] 0.931 -0.749 -0.47 -0.287 -0.216 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "1" "2" "3"
##   .. ..$ : chr [1:2] "Avg.Inc" "Total.Theft"
##  $ totss       : num 912
##  $ withinss    : num [1:3] 123.4 122.8 66.1
##  $ tot.withinss: num 312
##  $ betweenss   : num 600
##  $ size        : int [1:3] 197 220 40
##  $ iter        : int 2
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

Based on the model definition outputs, the resulting tot.withinss for the 3 years:

• 2018: 374

• 2019: 348

• 2020: 312

The resulting metrics are in line with the expectations set by the assault results. The overall range gives confidence the model comparisons are valid.

For 2018, the resulting cluster definitions for crime of theft:

• Cluster 1: Low income, low theft count

• Cluster 2: High theft count

• Cluster 3: High income, low theft count

For 2019, the resulting cluster definitions for crime of theft:

• Cluster 1: High theft count

• Cluster 2: Low income, low theft count

• Cluster 3: High income, low theft count

For 2020, the resulting cluster definitions for crime of theft:

• Cluster 1: High income, low theft count

• Cluster 2: Low income, low theft count

• Cluster 3: High theft count

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

## 825x669 terrain map image from Stamen Maps. 
## See ?ggmap to plot it.

As expected, the maps for each year are quite similar in location of higher theft crime counts. With the decreased number of clusters, the higher crime count cluster appears to be comprised of the equivalent of the medium and high crime count clusters from the total crime section. Also, the locations of the higher crime counts for assault and theft do appear similar. This assessment does not attempt to indicate one crime leads to another, but simply note the volume of assault and theft crimes do occur in similar locations in Cary.

References

• https://uc-r.github.io/kmeans_clustering

• https://www.tidymodels.org/learn/statistics/k-means/

• https://map-rfun.library.duke.edu/01_ggmap_georef.html

5 Principal Components Analysis

6 Geographic Weighted Regression

7 Conclusion

8 Appendices

8.1 References

National Incident-Based Reporting System (NIBRS) codes are described in the 2019 update: (https://www.fbi.gov/file-repository/ucr/ucr-2019-1-nibrs-technical-specification.pdf/view)

Vignettes for the use of sf package in R by Edzer Pebesma: (https://r-spatial.github.io/sf/articles/)

[Crime Prediction & Monitoring Framework Based on Spatial Analysis] (https://www.sciencedirect.com/science/article/pii/S187705091830807X)

Kounadi, Ristea, et al.(2020). A systematic review on spatial crime forecasting. Crime Science 9:7

8.2 Disclaimer

8.3 Code

We summarize all the R code used in this project in this appendix for ease of reading.

library(knitr)
# ---------------------------------
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)


# Your libraries go here

library(tidyverse) 
library(ggplot2) 
library(ggrepel)
library(kableExtra) 
library(dplyr)

library(caret)   # Model Framework
library(skimr)    # Used for EDA
library(klaR)     # Implemented KNN and Naive Bayes models, etc

# PLEASE ADD YOUR R LIBRARIES BELOW
# ------------------------------
library(janitor)  # used for tabyl function
library(sf)     # used for spatial data load and calculations
library(tmap)   # used for plotting maps


library(tidymodels)
library(cluster)
library(factoextra)
library(ggmap)
library(gridExtra)
# Alex:
# This code chunk is only to be used for conditionally disabling certain
# very slow model code as I merge different sections and need to frequent recompile
runSlowChunks = F


# Load the raw data for the town boundary
town_raw = sf::st_read("cary_boundary.geojson", quiet = TRUE)

# Load the raw data for the census
acs_raw = sf::st_read("acs2018_5yr_B19013_14000US37183053513.geojson", quiet = TRUE)

# Load the subdivisions
div_raw = sf::st_read("httpmapstownofcary0.geojson", quiet = TRUE)



cpd_raw = sf::st_read("cpd-incidents.geojson", quiet = TRUE)
dataset_sizes = data.frame( town_limits = as.character( dim(town_raw) ) ,
                            acs_income  = as.character( dim(acs_raw) ) ,
                            subdivisions = as.character( dim(div_raw) ),
                            crimes       = as.character( dim(cpd_raw) ) )

dataset_sizes = rbind( dataset_sizes, c("Town of Cary", "Censusreporter.org", "Town of Cary", "Town of Cary Police Department"))
row.names(dataset_sizes) = c("number rows", "number columns", "source")


dataset_sizes %>% kable( caption = "Raw data file dimensions" ) %>% kable_styling(bootstrap_options = c("striped", "hover"))
#
# Use the North Carolina NAD83 Projected CRS
# --------------------------------------------------
town = st_transform(town_raw, 2264)

# Remove any observations (subdivisions) where geometry data is missing.
#
div = div_raw %>% filter(  !is.na(st_dimension(div_raw) ) )

print(paste0("raw subdivisions data: ",  nrow(div_raw), " vs. ", "tidied subdivisions data: ", nrow(div)))

# Add an integer identifier to the subdivisions dataframe as there is no unique key (except geometry)
div = tibble::rowid_to_column(div, "div_id")

# Change the CRS coordinate system from WGS84 to NAD83/North Carolina
# in order to allow distance measurement, etc.
div = st_transform(div, 2264)


# Change the coordinate system for the ACS survey data to EPSG 2264 North Carolina ft-US
acs = st_transform(acs_raw, 2264)

acs = acs %>% rename(tract_name = name ) %>% 
  rename( income = B19013001 ) %>%  dplyr::select( -one_of(c("B19013001..Error") ) )

head(acs, 5) %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
summary(acs$income)
#  Calculation all the intersections between census tracts and
#  residential subdivisions
#
div_acs_partition = st_intersection( div, acs)

dim(div_acs_partition)
df_div_acs_areas = data.frame( area = as.numeric( st_area(div_acs_partition) ), 
                               tract_name = div_acs_partition$tract_name ,
                               name = div_acs_partition$name,
                               div_id = div_acs_partition$div_id, 
                               income = div_acs_partition$income ) 
df_div_acs_areas %>% group_by(div_id) %>%
                   summarize(avg_income = weighted.mean(income, area) ,
                             count = n(),
                             max_id = max(div_id) 
                             ) -> div_avg_income

div %>% dplyr::left_join( div_avg_income, by = 'div_id') %>% 
  dplyr::select( div_id, category, name, 
                  shape_stlength, shape_starea, approvedlots, description, geo_point_2d, 
                 avg_income, geo_point_2d) -> div

head(div, n = 3 ) %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover") )


# Only use the years 2018-2020 and drop those with missing or invalid geometry information.
#
cpd2018_2020 = cpd_raw %>% filter( !is.na(st_geometry(cpd_raw)))  %>%
                    filter(!is.na(st_dimension(cpd_raw))) %>%
                     filter(  year == "2018" | year == "2019" | year == "2020" ) # %>%  slice(1:10)

crimes = cpd2018_2020 %>% dplyr::select(incident_number, newdate, year, 
                                        ucr , crime_type, crime_category, offensecategory, 
                                        lon, lat, geocode, district)

# Change the coordinate system to NAD83/North Carolina ft US
crimes = st_transform(crimes, 2264)

dim(crimes)

head(crimes[c(1,3000,4000,6000), c("incident_number", "newdate", "ucr", "crime_type", "crime_category", 
                                   "geocode", "district")], n=4) %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
crimesdiv = st_join(crimes, div, join = st_within , left = FALSE)

dim(crimesdiv)

crimesdiv %>% group_by(year) %>% summarize(count=n())

sf::st_write(crimesdiv, "crimesdiv.geojson", append = FALSE, delete_dsn = TRUE, quiet = TRUE)
sf::st_write(div, "div.geojson", append = FALSE, delete_dsn = TRUE, quiet = TRUE)

crimesdiv = st_read("crimesdiv.geojson")

div = st_read("div.geojson")

crimesdiv = st_read("crimesdiv.geojson")

df.2018 <- crimesdiv %>% filter(year == "2018")
df.2019 <- crimesdiv %>% filter(year == "2019")

div_id_uniq <- unique(df.2018$div_id)
length(div_id_uniq) # for 2018, 485 unique div id
df.2018.1 <- df.2018 %>%
        group_by(div_id) %>%
        summarise(avg_inc = mean(avg_income), n = n())
sum(df.2018.1$n)
# n() gave the count, as each row is a "crime" ... but, if there are repeating incident numbers, let's try to pare down by unique incident number
#Year 2018
df.2018.2 <- df.2018 %>% as.data.frame %>%
        group_by(div_id) %>%
        summarise(avg_inc = mean(avg_income), n = length(unique(incident_number)))

df.2018.2$div_id <- as.character(df.2018.2$div_id)

# turn div_id into row names
df.2018.2 <- df.2018.2 %>%
         column_to_rownames('div_id')

#Year 2019
df.2019.2 <- df.2019 %>% as.data.frame %>%
        group_by(div_id) %>%
        summarise(avg_inc = mean(avg_income), n = length(unique(incident_number)))

df.2019.2$div_id <- as.character(df.2019.2$div_id)

# turn div_id into row names
df.2019.2 <- df.2019.2 %>%
         column_to_rownames('div_id')
sum(df.2018.2$n)
distance <- get_dist(df.2018.2)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
# Scale the data so the kmeans algorithm doesn't give too much weight to outliers
df.2018.2.prescale <- df.2018.2
df.2018.2 <- scale(df.2018.2)

df.2019.2.prescale <- df.2019.2
df.2019.2 <- scale(df.2019.2)
k2 <- kmeans(df.2018.2, centers=2, nstart=25)
str(k2)
fviz_cluster(k2, data=df.2018.2)
k3 <- kmeans(df.2018.2, centers = 3, nstart = 25)
k4 <- kmeans(df.2018.2, centers = 4, nstart = 25)
k5 <- kmeans(df.2018.2, centers = 5, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = df.2018.2) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = df.2018.2) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = df.2018.2) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = df.2018.2) + ggtitle("k = 5")

grid.arrange(p1, p2, p3, p4, nrow = 2)
set.seed(123)
fviz_nbclust(df.2018.2, kmeans, method = "wss")
set.seed(123)
fviz_nbclust(df.2018.2, kmeans, method = "silhouette")
set.seed(123)
gap_stat <- clusGap(df.2018.2, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
# Extracting results going with 5
set.seed(123)
final.18 <- kmeans(df.2018.2, 5, nstart = 25, algorithm="Hartigan-Wong")
final.19 <- kmeans(df.2019.2, 5, nstart = 25, algorithm="Hartigan-Wong")
# Extracting results going with 4
set.seed(123)
final.18_pam <- pam(df.2018.2, 4, metric = "euclidean", stand = FALSE)
# Basic plot
fviz_cluster(final.18_pam, data = df.2018.2)
clusplot(df.2018.2, final.18$cluster, color=T, shade=T, labels=1, lines=0)
# Basic plot
fviz_cluster(final.18, data = df.2018.2)
str(final.18)
df.2018.2.prescale <- as.data.frame(df.2018.2.prescale)
df.2018.2.prescale %>%
  mutate(Cluster = final.18$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")
df.2019.2.prescale <- as.data.frame(df.2019.2.prescale)
df.2019.2.prescale %>%
  mutate(Cluster = final.19$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")
dfX.18 <- data.frame('div_id' = names(final.18$cluster), 'cluster' = final.18$cluster)
rownames(dfX.18) <- NULL

dfX.19 <- data.frame('div_id' = names(final.19$cluster), 'cluster' = final.19$cluster)
rownames(dfX.19) <- NULL
# Join on div_id to add cluster by each subdivision

df.2018$div_id <- as.character(df.2018$div_id)
df.2018 <- left_join(df.2018, dfX.18, by='div_id')

df.2019$div_id <- as.character(df.2019$div_id)
df.2019 <- left_join(df.2019, dfX.19, by='div_id')
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.2018, aes(x=lon, y=lat, color=factor(cluster, labels = c("Low income, low crime count", "Medium crime count", "High crime count", "High income, low crime count", "Medium income, low crime count"))), size=0.5) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Total Crimes Clustered by Subdivision (2018)")

(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.2019, aes(x=lon, y=lat, color=factor(cluster, labels = c("Low income, low crime count", "High crime count", "Medium income, low crime count", "Medium crime count", "High income, low crime count"))), size=0.5) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Total Crimes Clustered by Subdivision (2019)")

df.18 <- crimesdiv %>% filter(year == "2018")
df.19 <- crimesdiv %>% filter(year == "2019")
df.20 <- crimesdiv %>% filter(year == "2020")
# Create wide DF to then capture the totals for theft and assault
# Year: 2018
df.18.wide <- spread(df.18, offensecategory, ucr)
colnames(df.18.wide)[which(names(df.18.wide) == "Aggravated Assault")] <- "Ag.As"
colnames(df.18.wide)[which(names(df.18.wide) == "Simple Assault")] <- "Sim.As"
colnames(df.18.wide)[which(names(df.18.wide) == "Stolen Property")] <- "Sto.Prop"
colnames(df.18.wide)[which(names(df.18.wide) == "Larceny")] <- "Larc"
colnames(df.18.wide)[which(names(df.18.wide) == "Motor Vehicle Theft")] <- "Mo.Ve.Th"
df.18.wide$Ag.As <- ifelse(is.na(df.18.wide$Ag.As), 0, 1)
df.18.wide$Sim.As <- ifelse(is.na(df.18.wide$Sim.As), 0, 1)
df.18.wide$Sto.Prop <- ifelse(is.na(df.18.wide$Sto.Prop), 0, 1)
df.18.wide$Larc <- ifelse(is.na(df.18.wide$Larc), 0, 1)
df.18.wide$Mo.Ve.Th <- ifelse(is.na(df.18.wide$Mo.Ve.Th), 0, 1)

df.19.wide <- spread(df.19, offensecategory, ucr)
colnames(df.19.wide)[which(names(df.19.wide) == "Aggravated Assault")] <- "Ag.As"
colnames(df.19.wide)[which(names(df.19.wide) == "Simple Assault")] <- "Sim.As"
colnames(df.19.wide)[which(names(df.19.wide) == "Stolen Property")] <- "Sto.Prop"
colnames(df.19.wide)[which(names(df.19.wide) == "Larceny")] <- "Larc"
colnames(df.19.wide)[which(names(df.19.wide) == "Motor Vehicle Theft")] <- "Mo.Ve.Th"
df.19.wide$Ag.As <- ifelse(is.na(df.19.wide$Ag.As), 0, 1)
df.19.wide$Sim.As <- ifelse(is.na(df.19.wide$Sim.As), 0, 1)
df.19.wide$Sto.Prop <- ifelse(is.na(df.19.wide$Sto.Prop), 0, 1)
df.19.wide$Larc <- ifelse(is.na(df.19.wide$Larc), 0, 1)
df.19.wide$Mo.Ve.Th <- ifelse(is.na(df.19.wide$Mo.Ve.Th), 0, 1)

df.20.wide <- spread(df.20, offensecategory, ucr)
colnames(df.20.wide)[which(names(df.20.wide) == "Aggravated Assault")] <- "Ag.As"
colnames(df.20.wide)[which(names(df.20.wide) == "Simple Assault")] <- "Sim.As"
colnames(df.20.wide)[which(names(df.20.wide) == "Stolen Property")] <- "Sto.Prop"
colnames(df.20.wide)[which(names(df.20.wide) == "Larceny")] <- "Larc"
colnames(df.20.wide)[which(names(df.20.wide) == "Motor Vehicle Theft")] <- "Mo.Ve.Th"
df.20.wide$Ag.As <- ifelse(is.na(df.20.wide$Ag.As), 0, 1)
df.20.wide$Sim.As <- ifelse(is.na(df.20.wide$Sim.As), 0, 1)
df.20.wide$Sto.Prop <- ifelse(is.na(df.20.wide$Sto.Prop), 0, 1)
df.20.wide$Larc <- ifelse(is.na(df.20.wide$Larc), 0, 1)
df.20.wide$Mo.Ve.Th <- ifelse(is.na(df.20.wide$Mo.Ve.Th), 0, 1)
  
df.18.v2 <- df.18.wide %>% as.data.frame %>%
        group_by(div_id) %>%
        summarise(Avg.Inc = mean(avg_income), 
                  Ag.Assault=sum(Ag.As),
                  Sim.Assault=sum(Sim.As), 
                  Sto.Prop=sum(Sto.Prop),
                  Larc=sum(Larc),
                  Mo.Ve.Th=sum(Mo.Ve.Th))

df.18.v2$div_id <- as.character(df.18.v2$div_id)

# turn div_id into row names
df.18.v2 <- df.18.v2 %>%
         column_to_rownames('div_id')

df.18.v2$Total.Assault <- df.18.v2$Ag.Assault + df.18.v2$Sim.Assault
df.18.v2$Total.Theft <- df.18.v2$Larc + df.18.v2$Sto.Prop + df.18.v2$Mo.Ve.Th

df.18.v2.TotAs <- df.18.v2[,c("Avg.Inc", "Total.Assault")]
df.18.v2.TotTh <- df.18.v2[,c("Avg.Inc", "Total.Theft")]

dim(df.18.v2.TotAs)
df.19.v2 <- df.19.wide %>% as.data.frame %>%
        group_by(div_id) %>%
        summarise(Avg.Inc = mean(avg_income), 
                  Ag.Assault=sum(Ag.As),
                  Sim.Assault=sum(Sim.As), 
                  Sto.Prop=sum(Sto.Prop),
                  Larc=sum(Larc),
                  Mo.Ve.Th=sum(Mo.Ve.Th))

df.19.v2$div_id <- as.character(df.19.v2$div_id)

# turn div_id into row names
df.19.v2 <- df.19.v2 %>%
         column_to_rownames('div_id')

df.19.v2$Total.Assault <- df.19.v2$Ag.Assault + df.19.v2$Sim.Assault
df.19.v2$Total.Theft <- df.19.v2$Larc + df.19.v2$Sto.Prop + df.19.v2$Mo.Ve.Th

df.19.v2.TotAs <- df.19.v2[,c("Avg.Inc", "Total.Assault")]
df.19.v2.TotTh <- df.19.v2[,c("Avg.Inc", "Total.Theft")]

dim(df.19.v2.TotAs)
df.20.v2 <- df.20.wide %>% as.data.frame %>%
        group_by(div_id) %>%
        summarise(Avg.Inc = mean(avg_income), 
                  Ag.Assault=sum(Ag.As),
                  Sim.Assault=sum(Sim.As), 
                  Sto.Prop=sum(Sto.Prop),
                  Larc=sum(Larc),
                  Mo.Ve.Th=sum(Mo.Ve.Th))

df.20.v2$div_id <- as.character(df.20.v2$div_id)

# turn div_id into row names
df.20.v2 <- df.20.v2 %>%
         column_to_rownames('div_id')

df.20.v2$Total.Assault <- df.20.v2$Ag.Assault + df.20.v2$Sim.Assault
df.20.v2$Total.Theft <- df.20.v2$Larc + df.20.v2$Sto.Prop + df.20.v2$Mo.Ve.Th

df.20.v2.TotAs <- df.20.v2[,c("Avg.Inc", "Total.Assault")]
df.20.v2.TotTh <- df.20.v2[,c("Avg.Inc", "Total.Theft")]

dim(df.20.v2.TotAs)
# Scale the data
# Save the data in prescale values for later use

df.18.v2.TotAs.Prescale <- df.18.v2.TotAs
df.18.v2.TotTh.Prescale <- df.18.v2.TotTh

df.18.v2.TotAs <- scale(df.18.v2.TotAs)
df.18.v2.TotTh <- scale(df.18.v2.TotTh)

df.19.v2.TotAs.Prescale <- df.19.v2.TotAs
df.19.v2.TotTh.Prescale <- df.19.v2.TotTh

df.19.v2.TotAs <- scale(df.19.v2.TotAs)
df.19.v2.TotTh <- scale(df.19.v2.TotTh)

df.20.v2.TotAs.Prescale <- df.20.v2.TotAs
df.20.v2.TotTh.Prescale <- df.20.v2.TotTh

df.20.v2.TotAs <- scale(df.20.v2.TotAs)
df.20.v2.TotTh <- scale(df.20.v2.TotTh)
# Total Assault
set.seed(123)

kclusts <- 
  tibble(k = 1:9) %>%
  mutate(
    kclust = map(k, ~kmeans(df.18.v2.TotAs, .x)),
    tidied = map(kclust, tidy),
    glanced = map(kclust, glance),
    augmented = map(kclust, augment, df.18.v2.TotAs)
  )

clusters <- 
  kclusts %>%
  unnest(cols = c(tidied))

assignments <- 
  kclusts %>% 
  unnest(cols = c(augmented))

clusterings <- 
  kclusts %>%
  unnest(cols = c(glanced))

p1 <- 
  ggplot(assignments, aes(x = Avg.Inc, y = Total.Assault)) +
  geom_point(aes(color = .cluster), alpha = 0.8) + 
  facet_wrap(~ k)

p2 <- p1 + geom_point(data = clusters, size = 10, shape = "x")
p2
ggplot(clusterings, aes(k, tot.withinss)) +
  geom_line() +
  geom_point()

df.18.v2.TotAs.final <- kmeans(df.18.v2.TotAs, centers=3, nstart=25)
str(df.18.v2.TotAs.final)

df.19.v2.TotAs.final <- kmeans(df.19.v2.TotAs, centers=3, nstart=25)
str(df.19.v2.TotAs.final)

df.20.v2.TotAs.final <- kmeans(df.20.v2.TotAs, centers=3, nstart=25)
str(df.20.v2.TotAs.final)

# Year: 2018
df.18.v2.TotAs.Prescale <- as.data.frame(df.18.v2.TotAs.Prescale)

df.18.v2.TotAs.Prescale %>%
  mutate(Cluster = df.18.v2.TotAs.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.18 <- data.frame('div_id' = names(df.18.v2.TotAs.final$cluster), 'cluster'=df.18.v2.TotAs.final$cluster)
rownames(dfX.18) <- NULL
# Year: 2019
df.19.v2.TotAs.Prescale <- as.data.frame(df.19.v2.TotAs.Prescale)

df.19.v2.TotAs.Prescale %>%
  mutate(Cluster = df.19.v2.TotAs.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.19 <- data.frame('div_id' = names(df.19.v2.TotAs.final$cluster), 'cluster'=df.19.v2.TotAs.final$cluster)
rownames(dfX.19) <- NULL
# Year: 2020
df.20.v2.TotAs.Prescale <- as.data.frame(df.20.v2.TotAs.Prescale)

df.20.v2.TotAs.Prescale %>%
  mutate(Cluster = df.20.v2.TotAs.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.20 <- data.frame('div_id' = names(df.20.v2.TotAs.final$cluster), 'cluster'=df.20.v2.TotAs.final$cluster)
rownames(dfX.20) <- NULL
# Join on div_id to add cluster by each subdivision
# Year: 2018
df.18.assault <- df.18 %>% filter(offensecategory == "Simple Assault" | offensecategory == "Aggravated Assault")

df.18.assault$div_id <- as.character(df.18.assault$div_id)
df.18.assault <- left_join(df.18.assault, dfX.18, by='div_id')

# Year: 2019
df.19.assault <- df.19 %>% filter(offensecategory == "Simple Assault" | offensecategory == "Aggravated Assault")

df.19.assault$div_id <- as.character(df.19.assault$div_id)
df.19.assault <- left_join(df.19.assault, dfX.19, by='div_id')

# Year: 2020
df.20.assault <- df.20 %>% filter(offensecategory == "Simple Assault" | offensecategory == "Aggravated Assault")

df.20.assault$div_id <- as.character(df.20.assault$div_id)
df.20.assault <- left_join(df.20.assault, dfX.20, by='div_id')
# Year: 2018
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.18.assault, aes(x=lon, y=lat, color=factor(cluster, labels = c("High Income, Low Assault Count", "High Assault Count", "Low Income, Low Assault Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Assaults Clustered by Subdivision (2018)")

# Year: 2019
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.19.assault, aes(x=lon, y=lat, color=factor(cluster, labels = c("High Income, Low Assault Count", "High Assault Count", "Low Income, Low Assault Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Assaults Clustered by Subdivision (2019)")

# Year: 2020
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.20.assault, aes(x=lon, y=lat, color=factor(cluster, labels = c("High Income, Low Assault Count", "High Assault Count", "Low Income, Low Assault Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Assaults Clustered by Subdivision (2020)")
# Total Theft
set.seed(123)

kclusts <- 
  tibble(k = 1:9) %>%
  mutate(
    kclust = map(k, ~kmeans(df.18.v2.TotTh, .x)),
    tidied = map(kclust, tidy),
    glanced = map(kclust, glance),
    augmented = map(kclust, augment, df.18.v2.TotTh)
  )

clusters <- 
  kclusts %>%
  unnest(cols = c(tidied))

assignments <- 
  kclusts %>% 
  unnest(cols = c(augmented))

clusterings <- 
  kclusts %>%
  unnest(cols = c(glanced))

p1 <- 
  ggplot(assignments, aes(x = Avg.Inc, y = Total.Theft)) +
  geom_point(aes(color = .cluster), alpha = 0.8) + 
  facet_wrap(~ k)

p2 <- p1 + geom_point(data = clusters, size = 10, shape = "x")
p2
ggplot(clusterings, aes(k, tot.withinss)) +
  geom_line() +
  geom_point()

df.18.v2.TotTh.final <- kmeans(df.18.v2.TotTh, centers=3, nstart=25)
str(df.18.v2.TotTh.final)

df.19.v2.TotTh.final <- kmeans(df.19.v2.TotTh, centers=3, nstart=25)
str(df.19.v2.TotTh.final)

df.20.v2.TotTh.final <- kmeans(df.20.v2.TotTh, centers=3, nstart=25)
str(df.20.v2.TotTh.final)

# Year: 2018
df.18.v2.TotTh.Prescale <- as.data.frame(df.18.v2.TotTh.Prescale)

df.18.v2.TotTh.Prescale %>%
  mutate(Cluster = df.18.v2.TotTh.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.18 <- data.frame('div_id' = names(df.18.v2.TotTh.final$cluster), 'cluster'=df.18.v2.TotTh.final$cluster)
rownames(dfX.18) <- NULL
# Year: 2019
df.19.v2.TotTh.Prescale <- as.data.frame(df.19.v2.TotTh.Prescale)

df.19.v2.TotTh.Prescale %>%
  mutate(Cluster = df.19.v2.TotTh.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.19 <- data.frame('div_id' = names(df.19.v2.TotTh.final$cluster), 'cluster'=df.19.v2.TotTh.final$cluster)
rownames(dfX.19) <- NULL
# Year: 2020
df.20.v2.TotTh.Prescale <- as.data.frame(df.20.v2.TotTh.Prescale)

df.20.v2.TotTh.Prescale %>%
  mutate(Cluster = df.20.v2.TotTh.final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")

dfX.20 <- data.frame('div_id' = names(df.20.v2.TotTh.final$cluster), 'cluster'=df.20.v2.TotTh.final$cluster)
rownames(dfX.20) <- NULL
# Join on div_id to add cluster by each subdivision
# Year: 2018
df.18.theft <- df.18 %>% filter(offensecategory == "Larceny" | offensecategory == "Stolen Property" | offensecategory == "Motor Vehicle Theft")

df.18.theft$div_id <- as.character(df.18.theft$div_id)
df.18.theft <- left_join(df.18.theft, dfX.18, by='div_id')

# Year: 2019
df.19.theft <- df.19 %>% filter(offensecategory == "Larceny" | offensecategory == "Stolen Property" | offensecategory == "Motor Vehicle Theft")

df.19.theft$div_id <- as.character(df.19.theft$div_id)
df.19.theft <- left_join(df.19.theft, dfX.19, by='div_id')

# Year: 2020
df.20.theft <- df.20 %>% filter(offensecategory == "Larceny" | offensecategory == "Stolen Property" | offensecategory == "Motor Vehicle Theft")

df.20.theft$div_id <- as.character(df.20.theft$div_id)
df.20.theft <- left_join(df.20.theft, dfX.20, by='div_id')
# Total thefts
# Year: 2018
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.18.theft, aes(x=lon, y=lat, color=factor(cluster, labels = c("Low Income, Low Theft Count", "High Theft Count", "High Income, Low Theft Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Thefts Clustered by Subdivision (2018)")

# Year: 2019
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.19.theft, aes(x=lon, y=lat, color=factor(cluster, labels = c("High Theft Count", "Low Income, Low Theft Count", "High Income, Low Theft Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Thefts Clustered by Subdivision (2019)")

# Year: 2020
(map <- get_map(c(left = -78.95, bottom = 35.64, right = -78.72, top = 35.87)))
ggmap(map) +
  geom_point(data=df.20.theft, aes(x=lon, y=lat, color=factor(cluster, labels = c("High Income, Low Theft Count", "Low Income, Low Theft Count", "High Theft Count"))), size=1) +
  xlab(expression(paste("Longitude (", degree,"W)"))) +
  ylab(expression(paste("Latitude (", degree,"N)"))) +
  labs(color = "Clusters") +
  ggtitle("Thefts Clustered by Subdivision (2020)")